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Abstract 

Lekner [J. Lekner, Mol. Simul. 20 (1998) ] and Sperb's [R. Sperb, Mol. Simul. 13, (1994)] 
work on the evaluation of Coulomb energy and forces under periodic boundary conditions is gen- 
eralized that makes it possible to use a triclinic unit cell in simulations in 3D rather than just an 
orthorhombic cell. The expressions obtained are in a similar form as previously obtained by Lekner 
and Sperb for the especial case of orthorhombic cell. 
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I. INTRODUCTION 



Molecular dynamics and Monte Carlo simulations of electrically charged point particles 
are indispensable in condensed matter physics. Generally, in simulations, periodic boundary 
conditions (pbc) are imposed to avoid unwanted effects of boundaries^. In presence of pbc, it 
becomes necessary to include the effect of image charges while calculating interaction energy 
and forces. The treatment of image charges is rather trivial for short range forces; one simply 
subdivides the simulation cell into several smaller cells, such that each of these cells is bigger 
than the range of interaction. The simulation time obviously scales as N for such short range 
potentials. However, if the interaction is long range, such as the Coulomb interaction, it 
becomes impossible to take into account all of image charges without taking recourse to 
some analytical technique. One of the techniques applied extensively in treating long range 
Coulomb interaction is the Ewald method^. Usually, when working with several thousand 
charges, the iV 2 cost of computing energy of a system carrying N charges overwhelms even 
supercomputers. To deal with such cases, one usually works with a variant of the Ewald 
method known as PPPM 3 . Another popular method is Fast Multipole method^. Recently 
another method known as the MMM^ method was shown to be faster and more accurate 
than the PPPM when a high accuracy is required, but has yet to attract the attention of 
researchers. Our aim in this paper is not to consider achieving linear scaling but rather to 
give a genuine alternative to the Ewald method for smaller systems. 

There are two main alternatives to the Ewald method. The first one is the so called 
Lekner method^, and the second one is due to Sperb^. These methods are limited in their 
application in that they were derived only for an orthorhombic simulation cell. In this 
regard, these methods are not as versatile as the Ewald method that can be applied even 
for a triclinic cell. However, only recently, a method was proposed 8 that extends Sperb's 
method and makes it possible to employ it even for a triclinic cell. In this paper, our aim is to 
directly generalize Lekner's work on orthorhombic simulation cells to a triclinic cell. Simple 
expressions will be derived to this end. However, our way of approaching this problem will 
be different from that of Lekner. The end results will of course contain Lekner's work as a 
special case. 
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II. GENERALIZATION OF LEKNER METHOD 



We consider N charges qi contained in a triclinic simulation cell in 3-dimensional space. 
The index i runs over i — 1 to N. We assume that the system is periodic in all dimensions. 
These charges interact via the Coulomb potential. The electrostatic energy of N charges 
can be expressed as 

^totai = g Yl ^1j G ( r i ~ r i) + 2 ^G'seif + y [J2 qiVi ) ' ( 2 ' X ) 

where the position of charges in the simulation cell is denoted by tv We will obtain expres- 
sions for G(r) and GWf in 3D in this section. 

The interaction between a pair of charges, a separation r apart, goes as | t* | . Such pair 
wise interactions when added under pbc lead to a diverging series, if the simulation cell is 
not overall charge neutral. However, if one has a charge neutral system, then the sum leads 
to a conditionally convergent series. To give a well defined meaning to the series, one has to 
specify how the terms in series are to be grouped together. Usually, one assumes that the 
particles interact with a screened potential that goes as exp (— (3 \ r\) / \r\; finally the limit 
(3 — > is taken. This is equivalent to introducing artificial background charges, and also 
to taking sums over expanding cubes. However, this technique only leads to the intrinsic 
part of the energy. A dipole term has to be added if one wants the energy in the limit of 
expanding spherical shells^. The last term in Eq. (j2.1J) represents this dipole term. 

We introduce a slightly different way of including the j3 factor. Instead of working with 
the exponential functions, we will work with the modified Bessel function of the second 
kind. Of course there is not much difference between these two different ways as for the 
large arguments, the modified Bessel functions of the second kind decay exponentially as 
well. We start with the fact that the limit of Ki/ 2 (/3r) as (3 tends to zero is given by 

to^M^-L^, (2.2) 

that makes it possible to write 

i../^ Ql/2 K 1/2 {(3r) 



lim/3 1 ^ : ' . (2.3) 

In a triclinic cell, the position of a charge can be specified by xi, x 2 and X3, where < Xi < k 
for i — 1,2, 3. Here Zj denote the lengths of the sides of the triclinic basic cell. To obtain 



the interaction between a pair of charges one may assume that one of the charges is located 
at the origin and the other one at (xi,x 2 ,x 3 ). Due to the pbc, one has to also consider 
the interaction of the second charge with all of the periodic images of the charge at the 
origin. These periodic images are located at {hm, l 2 n, l 3 p) , where m, n and p are integers 
ranging over — oo to +00. The distance between the first charge and a periodic image at 
(hm, l 2 n, l 3 p) is given by 

r m,n, P = ( x i + hm) 2 + (x 2 + l 2 nf + (x 3 + l 3 pf 
+ 2(xi + hm) (x 2 + l 2 n) cos a 
+ 2 (x 2 + l 2 n) (x 3 + l 3 p) cos j3 

+ 2 (x 3 + l 3 p) {xi + hm) cos 7. (2.4) 
The function G can be expressed now as 

G Or) = lim G (r; (3) 

= hm Jiff* V K ^\^—\ (2.5) 

m,n,p ' m,n,p 

Now, we switch over to the cylindrical coordinates by defining 

r2 m,n,p = (%np + h™? + pl, p , (2.6) 



where 



and 



x n ,p — xi + (x 2 + l 2 n) cos a + (x 3 + l 3 p) cos 7, (2.7) 



Pn, P — ( x 2 + hnf sin 2 a + (x 3 + hp) 2 sin 2 7 

+ 2 (#2 + ^) (^3 + hp) ( cos — cos a cos 7) . (2.8) 



For later convenience we also define 

x n, P = hn cos a + l 3 P cos 7, 



n,p 



(l 2 n) 2 sin 2 a + (/3P) 2 sin 2 7 + 2 (l 2 n) (l 3 p) (cos (3 — cos a cos 7) (2.9) 
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and 



Po,o — [ x l sm2 a + x l sm2 7 + 2a; 2 x 3 (cos /3 — cos a cos 7)] ^ 2 
^0,0 = Xi + X2 cos ct + £3 cos 7. 

Now, as shown in the appendix, G can be written as 

G{v) = y lim E ex P 6 27rm T^) ^° ( V^+^W) 
= [/(r) + Q(r), 



(2.10) 



where 



and 



n,p=— 00 



Q w = £ E E ^0 (2^^) cos 



(2.12) 



(2-13) 



n,p m=l 

We note that Q has excellent convergence for e = Po,o/h > 0.1. However, for smaller e we 
will modify Q. But before doing that we would like to obtain an expression for U. This can 
be done following Sperb. For this we first express p n , p as follows 

1 2 



P 2 n, P = sin2 a 



x 2 + l 2 n + (x 3 + l 3 p) 



cos /3 — cos a cos ' 
sin 2 a 



+ (a* + hp)' 



sin 2 7 



cos (5 — cos a cos 7 



sin a 



sin 2 a + l 2 n) 2 + (x 3 + hp) 2 Q 2 ] , 



(2.14) 



where 



y P = x 2 + (x 3 + Zap) C, 



(2.15) 



and 



c = 



cos P — cos a cos 7 



sin 2 a 



n = 



y/l — cos 2 a — cos 2 /? — cos 2 7 + 2 cos a cos /3 cos 7 

sin 2 a 



We note that ( and Q are purely geometrical factors. Also, note the relation 

,2 



fi 2 + C 2 



• 2 ' 
sm a 



(2.16) 
(2.17) 

(2.18) 
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which will be useful later. In order to further recast the expression for U we use the identity 



K ^(y + nlf + x* 



n exp (— [3 \x\) 
= 1 ]i 

exp(- W ^+(¥«) 2 ) 
+ V . = ^cos (27rny), 2.19 



which implies 

OO 

lim V K 

/3->0 



N(y + niy + x 2 



lim 7r exp(-/?|a;|) _ 1 
/3^o Z /3 2 



z' z. 



(2.20) 



where we have defined 



L[x, y) — In [1 — 2 exp (— 27r cos (27n/) + exp (— 47r |x| 



(2.21) 



Thus, we obtain from Eq. (l2~T3l and (l2~T9l 

2vr 1 °° 

U(r) = — KmM(x 3 ,{3)-- L 



p=—oo 



k ' k 



hk £-0 

The M in Eq. (J2~2~2l stands for 

M ( X3 ^) = ^ exp sina|x 3 + / 3 p|] 



(2.22) 



p=— OO 



QU sin a 



cosh 



e(i-2^ 



(2.23) 



2 £sinh£ 

where £ = (/?f2Z 3 /2) sin a; and a simple geometric sum has been carried out. The limit (3 — ► 
can now be carried out 



limM(x 3 ,/3) 



fiZ 3 sin a 



fiZ 3 sin a 



I_ 2 N 

3 h 



2|?M +i 
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1-6^+f^ 



/3 2 f2Z 3 sin a 



(2.24) 



Using Eq. (j2.11|) . ()2.22|) . ()2.23|) and ()2.24|) we obtain the following expression for the energy 
G (r) = i E E A-o (2™^) cos (a™^) - I £ L 



rt,p m=l 

f2Z 3 sin a 7r 
hk 3 



p=— 00 



^3 + ^3P n £/p 



'3 V '3 



(2.25) 
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where we have dropped the constant factor 47r/(V/3 2 ), and V = hhh fisina denotes the 
volume of the basic simulation cell. This dropping of the constant factor is justified on 
account of charge neutrality: The overall contribution to the energy from this term would 
be Ep 

-JV-l N / -i \ , N 



An 



V 







i=i j=i+i \r / i=1 \y / 



(2.26) 



The result in Eq. (J2.25j) is the main result of this paper. The form of G as written 
in Eq. (|2.25jl has excellent convergence for the most part of the simulation cell. However, 
the convergence is bad for the case when e = Po,o/h <C 1. The problem lies with the series 
corresponding to n = and p = in Eq. (j2.13|) for Q. For this case the argument of the 
function Kq becomes very small and it takes a lot of summation terms over m to achieve 
convergence. However, now there is now a well defined way to fix this problem; we isolate 
the series corresponding to n = and p = 0, and rewrite it in terms of Polygamma and Zeta 
functions*^. For this we first define 



f(x,P,l) 




[2.27) 



m=l 



Using Eq. (|2.25J) and 1)2.27]) . we can write 

^ n.p rn 1 % ' 

/ 



X3 + hP r, Vp 



'/a 



+ 2tt 



F3 



+ f(x 0l o,Po,o,h) + Y i { 2] n {^y)- L 



X3 n yo 
I, ' I, 



[2.28) 



where a prime over the summation sign indicates that the term corresponding to n = and 
p = is not to be included in the summation. As written in Eq. (), the series has two 
convergence problems when p ,o tends to zero. Firstly, when p .o — the last term is not 
defined. And secondly, as mentioned earlier, the function / (x, p, I) does not converge fast 
for small p. Both of these problem may be fixed by taking the limit po,o — * and combining 



the appropriate terms as follows. Using Eq. (j2.14|) we note that po,o —* means yo — > and 
x 3 — > 0. Also, 

P« iP = sin2 « [(%> + ^2^) 2 + (x 3 + l 3 p) 2 fi 2 ] 

implies 

Po,o = sin2 « [Vo + x l^ 2 ] ■ 
So, the last term in Eq. ()2.28|) may be written as 



£s ^0 
V' h. 



1 J / sin 2 a [y 2 + z|g j 2 
*i I V 4/ 2 ' 2 



2 

2n2i 



I| ln /'[?/o 2 + ^ 
^1 



7 2 

t 2 



^3 n yp_ 
«2 ' ^2 



— n — 

2 In 



h sin a 



(2.29) 
(2.30) 



(2.31) 



Now, through a simple Taylor expansion it could be shown 7 that for small y and z we have 

In (y 2 + z 2 )-L [y, z] = La [y, z] - 2 In (2vr) , 



where 

La[y, z) =2nz + ^- (y 2 - z 2 ) + ^ (y« - 6yV + J) 
2% 6 

H (y 6 — lhy A z 2 + lhy 2 z A — z 6 ) + higher order terms. 

2835 

With the help of identity in Eq. ()2.32j) we can write S for small po,o as follows: 



(2.32) 



S = —La 
h 



xs yo 



L\ II2 sina 



This fixes the first problem. To fix the second problem in order to achieve a better conver- 
gence for small p, we re-express / (x, p, I) using the results of Ref.— as follows: 

1 



f(x,p,l) 



r\ + r\ + rl + 2r\T2 cos a + 2r<ir?, cos (3 + 2r^r\ cos 7) 

jV-1 



1/2 



1 



1 



1 



' <•■ 1 \ 1/ p 2 + (n + x) z a/ p 2 4 (/' — .r) 

27 ^(N + z) + if>(N-x) 
1 



1 00 



m=l 



I 

"!/2 , „ 2m 



7?) 



p 2m [C (2m + 1, AT + x) + C (2m + 1, JV - a;)] 



(2.33) 
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where ip and ( stand for digamma and Hurwitz Zeta functions respectively, and iV > 1 is 
chosen such that iV > (p + x) . However, for better convergence it is desirable that iV > 
(p+ 1). Thus, using this alternate form for function / from Eq. (|2.33|) we obtain a different 
representation of G which gives very fast convergence as p tends to zero. The important 
fact here is that the Coulomb singularity toward small \r\ has been isolated. 
Also, it is now a simple matter to obtain the self-energy using the fact that 



lim S 

po, 0,3:3— >0 



lim I —La 



^3 Q Vo 



An 



I2 sin a 



An 



l 2 sin a 



It then follows 



G, 



self 



lim 



G(r) 



^ ' p 2 + x 2 



A 1 00 / 

TV K a 2nm Pn * 



n,p m=l 



hP a (hp 



x cos ( 27rm— 



An 



li \l 2 sma 



+ 



27 
h 



(2.34) 



where 



Pn.r. 



sin 2 a [(l 3 p( + l 2 nf + (hp) 2 W}. 



Expressions in Eq. ()2.25|) and ()2.34|) are the generalization of Sperb's and Lekner's work 
from an orthorhombic cell to a a triclinic cell. In the special case when a orthorhombic cell 
is considered, Eq. (|2.25|) reduces to 



G 



ortho 



(x 2 + hnf + (x 3 + hpf 



n,p m=l 



00 



p=— 00 



1 — 2 exp 



2tt 



£3 + hp\ 



cos 



( 2nm^- 



%2 



x cos ( 2n— + exp I — — \x 3 + hp\ 



An 



+ - In ( x /x 2 2 + x% ) + / ( x u \jx 2 2 + x£, h 



2n 



\x 3 \ 



+ 2n 



.1. ■> 



hhh 



[2.35) 
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This result is in agreement with the ones derived by Sperb and Lekner. Finally, to summarize 
our results, the total energy of iV charges is given by Eq. (j2.1j) . The function G in Eq. (j2.1|) 
may be obtained by using 



G (r) = ± £ t *o (2™^) - (2™^) - 1 f L 

n,p m=l x ' x ' p= — 00 



3^3 + ^3P n % 
/ 2 ' k 



+ 



fiZa sin a n 



hh 



(2.36) 



where x„ iP , y p are defined in Eqs. ()2.8|) . ()2.7j) . and (|2.15j) . The is defined by Eq. 
()2.17|) and the function L is defined in Eq. (j2.21|) . For the case when e = Po,o/h <C 1, one 
should use the following expression 



G (r) = i V V K ( 2nm^A cos f 27rm^ 



*1 



3?3 + Z 3 P n 2/p 
Z 2 ' h 



2tt 



F3 



Z1Z2Z3 hh 



27 1 

+ / (x ,o, po,o, l i) + ~r + i-La 
h h 



47T 



^2.37) 



/1 \l 2 sina / 

where and / is is defined in Eq. f!2.33|) and La is defined in Eq . (I2.32J) . The expression for 
force may be obtained by using F (r) = -VG (r). Finally, the self-energy required in Eq. 
(12.1 J) is obtained by using 



G 



self 



n,p m=l 



P 



7rm- 



n,p 
IT 



X COS I 27T771 



hp n (hp 



47T 



V,,. 

<i \<2 sin a 



(2.38) 



III. RESULTS AND CONCLUSIONS 



We have obtained complete expressions for the Coulomb potential in 3D, including self- 
energies. The results were derived for the most general case of a triclinic cell in 3D. The 
formulas derived here are in the same form as derived earlier by Lekner— and Sperb^ for 
the case of an orthorhombic cell. With our results, it should now be possible to extend any 
written code for Lekner or Sperb summation to a triclinic cell with minimal effort. The 
results obtained in this paper reduces to the results of a recent paper— when all angles 



10 



pertaining to the unit cell are set to ir/2. The results also agree with the special cases 
considered by Lekner— and Sperb*^. 



APPENDIX A: REPRESENTATION FOR G(r; 0) 



The solution of 

{V 2 -p 2 )S d (r;P) = -CMr) 
in ci-dimensional space is given by 



S d (r;P) 



0' 



,K v {(5r) 



(Al) 



(A2) 



C,= I 



where Cd is defined as 

2 d = l 

2tt d = 2, 

An u+1 /T (v) d>2. 

Here, T(u) stands for the Gamma function, and v = (d — 2) /2. For d = 3 we can write the 
solution in the Fourier space using Eq. (|Alj) : 

„i2n(k-i x+kn.p) 

(A3) 



S 3 (r; j3) = n9 / dkidk 



where we have written the vector r in cylindrical coordinates as (x, p) . Now, using Eq. ()2.5|) 
we can write 



c 



\3/2 



Kin (/Jr. 



/2 IP" m,n,p) 



,1/2 



(2*r' 



(27T) 



E e / « fe 



3 i27r[fei(a;„ iJ ,+Zim)+fep.p)] 



^ 1.2 



A; 2 + k l p + /? 2 /4tt 2 



(A4) 



Now, using the fact that 



m) 



(A5) 



we can replace the integral over k% with a summation: 
G(r;0) 



(7 I /■ i27r(mx n ,p+fc p .p n >p ) 

" EE / dk » 



m 2 + fc 2 + /? 2 /4tt 2 



S E 6X P O 2 ™^) K ° {V^W 2 Pn, P ) , (A6) 



7l,p 111 
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where 

Pn,p — \Pn,p\ 
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